Introducción

Veíamos en los modelos gSMR que hay una separación perfecta entre animales marcados y no marcados (recordad que los eventos de marcaje desconocido, cuando no se puede discriminar si los animales son marcados o no, se desechan). Sin embargo, un problema común es que para los animales marcados no se puedan leer las marcas, y se pierda por tanto la identificación. Es decir, aun sabiendo que están marcados, no se puede reconocer el animal individualmente. Para abordar este problema, Jiménez et al. (2019) desarrollaron el marcaje-reavistamiento espacial generalizado con identificación imperfecta (gSMR-ID). Al igual que rtSCR se basa en un muestreador Metropolis-Hastings (MH) que impone una restricción en cada iteración MCMC que consiste en que el total de eventos de individuos marcados debe ser igual a la suma de eventos de individuos marcados e identificados (conocidos) y de eventos de marcados y no identificados (cuya identidad es latente) y para los cuales solo observamos eventos de conteos nnid[j,k]. Pongo unas fotografías de cada tipología para aclarar conceptos.



Marcados identificados





Marcados no identificados





No marcados





Evento de marcaje desconocido (a desechar)



En este modelo, en vez de emplear simulaciones para generar los datos, vamos a usar datos reales de zorro (Vulpes vulpes) en el Parque Nacional de Cabañeros (centro de España). Para ejecutar este modelo, debeis descargar los datos en zip (DataZorro.zip), descomprimirlos y ponerlos en vuestro directorio de trabajo. A continuacion, ejecutad el modelo. Se trata de un caso extremo, donde ninguno de los zorros marcados era identificable, ya que el patrón de marcaje con bandas en marcas auriculares con bandas resultó ineficaz para el reconocimiento con fototrampeo.

En el proceso de marcaje se marcaron 11 animales. En el proceso de reavistamiento obtuvimos 28 fotos con estado de marcaje desconocido (a desechar); 29 eventos de marcados sin identificar; ningún de marcado identificado, y 396 eventos sin marcaje. Marcamos con GPS 4 individuos. Una cuestión a señalar es que los procesos de marcaje y de reavistamiento se solapaban. Por ello usamos una máscara binaria (0/1) con la fecha de marcaje, a fin de acotar las probabilidades de reavistamiento a las fechas en que los animales estaban efectivamente marcados. Ejecutamos el modelo usando R (R Core Team, 2022) y Nimble (De Valpine et al., 2017).

Preparamos los datos

# El número total de ocasiones de muestreo utilizados fue de 21 dias
# (06/3/2014-30/3/2014)

# Operatividad del proceso de marcado
Oper.mark<- data.matrix(read.table("Oper_trap.txt", header=FALSE))
dim(Oper.mark)
## [1] 31 21
Oper.mark.R<-apply(Oper.mark,1,sum)

# PROCESO DE MARCAJE
#====================
mark.ch <- secr::read.capthist("ymark.txt", "trapssecr.txt", noccasions =21)
## No errors found :-)
str(mark.ch)
##  'capthist' num [1:11, 1:21, 1:31] 0 0 0 0 0 0 0 0 0 0 ...
##  - attr(*, "dimnames")=List of 3
##   ..$ : chr [1:11] "VvC01" "VvC02" "VvC03" "VvC04" ...
##   ..$ : chr [1:21] "1" "2" "3" "4" ...
##   ..$ : chr [1:31] "1" "2" "3" "4" ...
##  - attr(*, "covariates")='data.frame':   0 obs. of  0 variables
##  - attr(*, "traps")=Classes 'traps' and 'data.frame':    31 obs. of  2 variables:
##   ..$ x: int [1:31] 384405 383724 382697 381054 380908 380831 384279 384082 382639 383014 ...
##   ..$ y: int [1:31] 4353660 4354856 4356469 4356155 4355523 4351881 4353536 4353486 4356457 4354875 ...
##   ..- attr(*, "detector")= chr "multi"
##   ..- attr(*, "spacex")= num 5
##   ..- attr(*, "spacey")= num 1
##   ..- attr(*, "spacing")= num 133
##  - attr(*, "session")= chr "1"
##  - attr(*, "inject.time")= num [1:11] 0 0 0 0 0 0 0 0 0 0 ...
y.mark3D<-as.array(mark.ch)
y.mark<-apply(y.mark3D,c(1,3),sum)
dim(y.mark)
## [1] 11 31
(n.marked<-dim(y.mark)[1]) # Número de marcados
## [1] 11
# Testamos si está correcto el archivo de operatividad del marcaje
y.mark3D.test<-aperm(y.mark3D, c(1,3,2))
sum(apply(y.mark3D.test,c(2,3),sum)); dim(y.mark3D.test)
## [1] 11
## [1] 11 31 21
sum(apply(y.mark3D.test,c(2,3),sum)* Oper.mark)
## [1] 11
which(apply(y.mark3D.test,c(2,3),sum)* 
        Oper.mark-apply(y.mark3D.test,c(2,3),sum)<0, arr.in=TRUE)
##      row col
# Trampas utilizadas
(J.mark<-dim(Oper.mark)[1])
## [1] 31
# Dias de marcaje
(K.mark<- dim(Oper.mark)[2])
## [1] 21
M<-750 # Aumentado de datos

y.mark.aug<-array(0,c(M,J.mark))
y.mark.aug[1:11,]<-y.mark

# Ploteamos la operatividad
x1<-as.matrix(Oper.mark)  ## Tenemos que convertir esto en matriz
image(1:K.mark,1:J.mark,t(x1), yaxt = "n", xlab="Ocasión", ylab="", 
      cex.lab=1.25, col=topo.colors(2))
mtext(side = 2, "Trampas en vivo", line = 2.5, cex=1.25)
axis(2, rev(seq(1, J.mark, by=2)))
box()

# Fechas de marcaje (para crear la máscara binaria)
(mark.date<-data.matrix(read.table("mark_date.txt", header=FALSE)))
##       V1 V2 V3 V4 V5 V6 V7 V8 V9 V10 V11 V12 V13 V14 V15 V16 V17 V18 V19 V20
##  [1,]  0  1  1  1  1  1  1  1  1   1   1   1   1   1   1   1   1   1   1   1
##  [2,]  0  1  1  1  1  1  1  1  1   1   1   1   1   1   1   1   1   1   1   1
##  [3,]  0  0  0  1  1  1  1  1  1   1   1   1   1   1   1   1   1   1   1   1
##  [4,]  0  0  0  1  1  1  1  1  1   1   1   1   1   1   1   1   1   1   1   1
##  [5,]  0  0  0  1  1  1  1  1  1   1   1   1   1   1   1   1   1   1   1   1
##  [6,]  0  0  0  0  0  1  1  1  1   1   1   1   1   1   1   1   1   1   1   1
##  [7,]  0  0  0  0  0  0  1  1  1   1   1   1   1   1   1   1   1   1   1   1
##  [8,]  0  0  0  0  0  0  1  1  1   1   1   1   1   1   1   1   1   1   1   1
##  [9,]  0  0  0  0  0  0  1  1  1   1   1   1   1   1   1   1   1   1   1   1
## [10,]  0  0  0  0  0  0  0  0  0   0   0   0   0   0   0   1   1   1   1   1
## [11,]  0  0  0  0  0  0  0  0  0   0   0   0   0   0   0   0   0   0   0   0
##       V21 V22 V23 V24 V25 V26 V27 V28 V29 V30 V31 V32 V33 V34 V35 V36 V37 V38
##  [1,]   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1
##  [2,]   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1
##  [3,]   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1
##  [4,]   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1
##  [5,]   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1
##  [6,]   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1
##  [7,]   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1
##  [8,]   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1
##  [9,]   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1
## [10,]   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1
## [11,]   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1
##       V39 V40 V41 V42 V43 V44 V45 V46 V47 V48 V49
##  [1,]   1   1   1   1   1   1   1   1   1   1   1
##  [2,]   1   1   1   1   1   1   1   1   1   1   1
##  [3,]   1   1   1   1   1   1   1   1   1   1   1
##  [4,]   1   1   1   1   1   1   1   1   1   1   1
##  [5,]   1   1   1   1   1   1   1   1   1   1   1
##  [6,]   1   1   1   1   1   1   1   1   1   1   1
##  [7,]   1   1   1   1   1   1   1   1   1   1   1
##  [8,]   1   1   1   1   1   1   1   1   1   1   1
##  [9,]   1   1   1   1   1   1   1   1   1   1   1
## [10,]   1   1   1   1   1   1   1   1   1   1   1
## [11,]   1   1   1   1   1   1   1   1   1   1   1
mark.date.aug<-array(0,c(M,49))
mark.date.aug[1:11,]<-mark.date


# PROCESO DE REAVISTAMIENTO
#===========================
# Operatividad de las cámaras-trampa
Oper.resight<-  data.matrix(read.table("Oper_cam.txt", header=FALSE))
(K <- dim(Oper.resight)[2])
## [1] 49
(J.resight <- dim(Oper.resight)[1])
## [1] 40
x2<-as.matrix(Oper.resight)  ## Tenemos que convertir esto en matriz
image(1:K,1:J.resight,t(x2), yaxt = "n", xlab="Ocasión", ylab="", 
      cex.lab=1.25, col=topo.colors(2))
mtext(side = 2, "Cámaras-trampa", line = 2.5, cex=1.25)
axis(2, rev(seq(1, J.resight, by=2)))
box()

# Marcados observados con identificación (ninguno)
datYknown<-array(0,c(n.marked,J.resight,K))
sum(datYknown)
## [1] 0
# Marcados observados sin identificación
nniDD<- secr::read.capthist("nnid.txt", "traplocs.txt", noccasions =49)
## No errors found :-)
nniDD<-aperm(nniDD,c(1,3,2))
nnid<-data.matrix(apply(nniDD,c(2,3),sum))
sum(nnid)  # marcados sin identificar
## [1] 29
# Para buscar los valores negativos (errores en Oper.resight) usamos:
which((nnid*Oper.resight)-nnid<0, arr.in=TRUE)
##      row col
# Escalamos y centramos coordenadas de las cámaras
x.resight<-data.matrix(read.table("traplocs.txt", header=FALSE))[,2:3]/1000
Xx<-mean(x.resight[,1])
Xy<-mean(x.resight[,2])
x.resight[,1]<-(x.resight[,1]-Xx)
x.resight[,2]<-(x.resight[,2]-Xy)

# Escalamos y centramos coordenadas de las trampas
livetraps<-as.matrix(secr::traps(mark.ch))
x.mark<-data.matrix(livetraps)/1000
x.mark[,1]<-(x.mark[,1]-Xx)
x.mark[,2]<-(x.mark[,2]-Xy)

# Eventos sin marcar
n<-secr::read.capthist("unmarked.txt", "traplocs.txt", noccasions =49)
## No errors found :-)
n<-data.matrix(apply(n,c(3,2),sum))
sum(n)
## [1] 396
# Para testar los valores negativos (errores en Oper.resight) usamos:
which(n*Oper.resight -n<0, arr.in=TRUE)
##      row col
# Gráfico de capturas
plot(x.resight,
     xlim=c(min(x.resight[,1])- 1,max(x.resight[,1])+ 1),
     ylim=c(min(x.resight[,2])- 1,max(x.resight[,2])+ 1),
     pch="+", 
     xlab="X",
     ylab="Y",
     asp=1)
points(x.mark, pch="+", col="red")
tot<-apply(nnid, 1,sum)
symbols(x.resight, circles=tot/10, inches=F,bg="#EEC90033", fg=NULL, add=T)
tot2<-apply(n,1,sum)
symbols(x.resight, circles=tot2/10, inches=F,bg="#228B2219", fg=NULL, add=T)


# DATOS DE TELEMETRÍA (4 IND)
#============================
locs<-read.table("PosZorro.txt", header=FALSE)[,1:2]/1000
locs[,1]<-locs[,1]-Xx
locs[,2]<-locs[,2]-Xy
locs1<-locs[1:31,]
locs2<-locs[32:461,]
locs3<-locs[462:527,]
locs4<-locs[528:546,]

points(locs1[,1], locs1[,2], pch=16, col="#FF00004C", cex=1)
points(locs2[,1], locs2[,2], pch=16, col="#0000FF4C", cex=1)

xy1<-cbind(locs1[,1], locs1[,2])
xy1<-data.matrix(xy1)
xysp1 <- SpatialPoints(xy1)
kudl1 <- kernelUD(xysp1, grid=100)
homerange1 <- getverticeshr(kudl1); homerange1@data$area*1E6
## [1] 126.5946
contour(getvolumeUD(kudl1)[1], level=c(50,95), add=TRUE, col=c('blue','blue'), lwd=2)

xy2<-cbind(locs2[,1], locs2[,2])
xy2<-data.matrix(xy2)
xysp2 <- SpatialPoints(xy2)
kudl2 <- kernelUD(xysp2, grid=100)
homerange2 <- getverticeshr(kudl2); homerange2@data$area*1E6
## [1] 706.2289
contour(getvolumeUD(kudl2)[1], level=c(50,95), add=TRUE, col=c('blue','blue'), lwd=2)

points(locs3[,1], locs3[,2], pch=16, col="#00FF004C", cex=1)
points(locs4[,1], locs4[,2], pch=16, col="forestgreen", cex=1)

locs<-rbind(locs1,locs2,locs3,locs4)
ind<-c(rep(1,nrow(locs1)),rep(2,nrow(locs2)),rep(3,nrow(locs3)),rep(4,nrow(locs4)))
nlocs<-nrow(locs)
n.collar<-4

# ESPACIO DE ESTADOS
buff<- 2
xl<-min(x.resight[,1])-buff
xu<-max(x.resight[,1])+buff
yl<-min(x.resight[,2])-buff
yu<-max(x.resight[,2])+buff
xlims=c(xl, xu)
ylims=c(yl, yu)
area<-(xu-xl)*(yu-yl)

Modelo en nimble

NimModel <- nimbleCode({

  lam0.resight ~ dunif(0,2) # Probabilidad basal de reavistamiento
  lam0.mark  ~ dunif(0,1)   # Probabilidad basal de captura
  sigma ~ dunif(0,2)        # parámetro de forma de la seminormal
  sigma2 <- sigma^2
  id.prob ~ dunif(0,1)      # probabilidad de identificación de marcados
  psi ~ dbeta(1,1)          # parámetro del aumentado de datos

  for(i in 1:M) {
    z[i] ~ dbern(psi)
    s[i,1] ~ dunif(xlim[1], xlim[2]) # límites x del espacio de estados
    s[i,2] ~ dunif(ylim[1], ylim[2]) # límites x del espacio de estados
    # Cálculos auxiliares para el marcado
    d.mark2[i,1:J.mark] <- (s[i,1]-x.mark[1:J.mark,1])^2 + 
                           (s[i,2]-x.mark[1:J.mark,2])^2
    lam.mark[i,1:J.mark] <- lam0.mark*
                            exp(-d.mark2[i,1:J.mark]/(2*sigma2))*z[i]
     # Cálculos auxiliares para el reavistamiento
    d.resight2[i,1:J.resight] <- (s[i,1]-x.resight[1:J.resight,1])^2 + 
                                 (s[i,2]-x.resight[1:J.resight,2])^2
    lam[i,1:J.resight] <- lam0.resight * 
                          exp(-d.resight2[i,1:J.resight]/(2*sigma2))*z[i]
    
    # Proceso de marcado
    for(j in 1:J.mark) {
      y.mark[i,j] ~ dbinom(lam.mark[i,j], Oper.mark.R[j]) 
    }
  }
  # Telemetría para los 4 animales con GPS
  for(r in 1:n.locs){
    locs[r,1]~dnorm(s[ind[r],1], 1/(sigma2))
    locs[r,2]~dnorm(s[ind[r],2], 1/(sigma2))
  }

  # Reavistamiento de los animales marcados (ninguno identificado)
  for(i in 1:n.marked)  {
    for(j in 1:J.resight) {
      for(k in 1:K){
        # Modelo para los historiales de captura completos de los marcados
        # considerando la fecha de marcaje
        y.full[i,j,k] ~ dpois(lam[i,j]*mark.date[i,k]*Oper.resight[j,k]) 
        # Modelo para la identificación de los marcados
        y.obs[i,j,k] ~ dbin(id.prob, y.full[i,j,k])
      }
    }
  }

  # Conteos de los animales no marcados
  for(j in 1:J.resight) {
    for(k in 1:K){
      # Consideramos la fecha de marcaje
      Lam[j,k] <- sum(lam[1:M,j]*(1-mark.date[1:M,k]))
      n[j,k] ~ dpois(Lam[j,k]*Oper.resight[j,k])
    }
  }


  N <- sum(z[1:M])
  D<-N/area
 })

Muestreador Metropolis-Hastings

## Muestreador MH para actualizar y[1:n.marked,j,k] de forma que 
## sume siempre nnid[j,k] + y.obs[i,j,k]
IDSampler <- nimbleFunction(
  contains = sampler_BASE,
  setup = function(model, mvSaved, target, control) {
    # Definimos componentes
    nnid<-control$nnid
    j<-control$j
    k<-control$k
    n.marked<-control$n.marked
    #nnid<-control$nnid
    calcNodes <- model$getDependencies(target)
  },

  run = function() {
    # conteos esperados de individuos por trampa
    lam.curr <- model$lam[1:n.marked,j]

    # Muestreo de y[1:n.marked,j] reasignando n[j] usando la
    # probabilidad condicional completa
    switch.probs <- lam.curr[1:n.marked]/sum(lam.curr[1:n.marked])

    # propone nueva identificaciones para nnid[j,k]
    y.latent.curr <- model$y.full[1:n.marked,j,k] - model$y.obs[1:n.marked,j,k]
    y.latent.prop <- rmulti(1, nnid, switch.probs[1:n.marked])
    model$y.full[1:n.marked,j,k] <<- model$y.obs[1:n.marked,j,k] + y.latent.prop

    # Modelo inicial logProb
    model_lp_initial <- model$getLogProb(calcNodes)

    # Modelo propuesto logProb
    model_lp_proposed <- model$calculate(calcNodes)

    # log-Metropolis-Hastings ratio
    log_MH_ratio <- (model_lp_proposed + dmulti(y.latent.curr, nnid, 
                                                switch.probs, log=TRUE)) -
                    (model_lp_initial + dmulti(y.latent.prop,  nnid, 
                                               switch.probs, log=TRUE))

    # Pasos Metropolis-Hastings
    accept <- decide(log_MH_ratio)
    if(accept) {
      copy(from = model, to = mvSaved, row = 1, nodes = calcNodes, 
           logProb = TRUE)
    } else {
      copy(from = mvSaved, to = model, row = 1, nodes = calcNodes, 
           logProb = TRUE)
    }
  },
  methods = list( reset = function () {} )
)

Constantes

constants      <- list(nnid=nnid,
                       J.resight=J.resight,
                       J.mark=J.mark,
                       M=M,
                       K=K,
                       xlim=xlims,
                       ylim=ylims,
                       n.marked=n.marked,
                       n.collar=4,
                       n.locs=nlocs,
                       ind=ind,
                       area=area)
str(constants)
## List of 12
##  $ nnid     : num [1:40, 1:49] 0 0 0 0 0 0 0 0 0 0 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:40] "1" "2" "3" "4" ...
##   .. ..$ : chr [1:49] "1" "2" "3" "4" ...
##  $ J.resight: int 40
##  $ J.mark   : int 31
##  $ M        : num 750
##  $ K        : int 49
##  $ xlim     : num [1:2] -10 10.2
##  $ ylim     : num [1:2] -5.44 6.05
##  $ n.marked : int 11
##  $ n.collar : num 4
##  $ n.locs   : int 546
##  $ ind      : num [1:546] 1 1 1 1 1 1 1 1 1 1 ...
##  $ area     : num 233

Datos

data     <-    list   (y.obs=datYknown,
                       y.mark=y.mark.aug,
                       n=n,
                       locs=locs,
                       x.resight=x.resight,
                       x.mark=x.mark,
                       mark.date=mark.date.aug,
                       Oper.resight=Oper.resight,
                       Oper.mark.R=Oper.mark.R)
                
str(data)
## List of 9
##  $ y.obs       : num [1:11, 1:40, 1:49] 0 0 0 0 0 0 0 0 0 0 ...
##  $ y.mark      : num [1:750, 1:31] 0 0 0 0 0 0 0 0 0 0 ...
##  $ n           : num [1:40, 1:49] 0 0 0 0 0 0 0 0 0 0 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:40] "1" "2" "3" "4" ...
##   .. ..$ : chr [1:49] "1" "2" "3" "4" ...
##  $ locs        :'data.frame':    546 obs. of  2 variables:
##   ..$ V1: num [1:546] 4.32 4.34 4.23 4.3 4.11 ...
##   ..$ V2: num [1:546] 0.907 1.09 0.84 0.895 0.979 ...
##  $ x.resight   : num [1:40, 1:2] -4.9 -4.06 -2.23 5.88 -7.15 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr [1:2] "V2" "V3"
##  $ x.mark      : num [1:31, 1:2] 7.28 6.6 5.57 3.93 3.78 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:31] "CO01" "CO02" "CO03" "CO04" ...
##   .. ..$ : chr [1:2] "x" "y"
##  $ mark.date   : num [1:750, 1:49] 0 0 0 0 0 0 0 0 0 0 ...
##  $ Oper.resight: int [1:40, 1:49] 1 0 1 1 1 1 1 1 1 1 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr [1:49] "V1" "V2" "V3" "V4" ...
##  $ Oper.mark.R : int [1:31] 8 8 3 2 2 6 4 3 5 3 ...

Inicios

# Marcados ID
s.start <- cbind(runif(M, xlims[1], xlims[2]), runif(M, ylims[1], ylims[2]))
d <- e2dist(s.start[1:n.marked,], x.resight)
lam0s<- 0.1
sigs <- 0.8
lam <- lam0s * exp( -(d^2)/(2 * sigs^2))
lam[lam=="NaN"]<-0

# Marcados no ID
yi <- array(0, c(n.marked, J.resight, K))
for (j in 1:J.resight) {
  for (k in 21:K) {  # uso a partir de 21 (todos marcados)
    if (nnid[j, k] > 0) {
      probs <- lam[ ,j]
      probs <- probs / sum(probs)
      latent.id <- sample(1:n.marked, nnid[j,k], prob = probs, replace = FALSE)
      yi[latent.id , j, k] <- 1
    }
  } # end of k
}   # end of j
datYknown<-0
yi <- yi + datYknown
apply(yi,1,sum)  # ejecutar el inicio hasta sum(yi)==sum(apply(yi,c(1,3),sum)*mark.date)
##  [1] 0 4 0 0 6 0 0 0 0 0 0
sum(yi)
## [1] 10
sum(apply(yi,c(1,3),sum)*mark.date)
## [1] 10
zst<-rep(1,M)
id.prob.s<-sum(datYknown)/(sum(datYknown)+sum(nnid))
n.collar<-4

inits   <-        list(z=zst,
                       s=s.start,
                       lam0.resight=0.1,
                       lam0.mark=0.002,
                       sigma=sigs,
                       id.prob=id.prob.s,
                       psi=0.2,
                       y.full=yi)
str(inits)
## List of 8
##  $ z           : num [1:750] 1 1 1 1 1 1 1 1 1 1 ...
##  $ s           : num [1:750, 1:2] -4.47 2.49 5.15 7.83 6.61 ...
##  $ lam0.resight: num 0.1
##  $ lam0.mark   : num 0.002
##  $ sigma       : num 0.8
##  $ id.prob     : num 0
##  $ psi         : num 0.2
##  $ y.full      : num [1:11, 1:40, 1:49] 0 0 0 0 0 0 0 0 0 0 ...

Parámetros a muestrear

params <- c('psi', 'lam0.resight', 'lam0.mark', 'sigma', 'N', 'D', 'id.prob')

Ejecutamos el modelo

start.time<-Sys.time()
Rmodel <- nimbleModel(code=NimModel,
                      constants=constants,
                      data=data,
                      inits=inits,
                      check=FALSE)
## Defining model
## Building model
## Setting data and initial values
## Running calculate on model
##   [Note] Any error reports that follow may simply reflect missing values in model variables.
## Checking model sizes and dimensions
Rmodel$initializeInfo()
##   [Note] All model variables are initialized.
Cmodel <- compileNimble(Rmodel)
## Compiling
##   [Note] This may take a minute.
##   [Note] Use 'showCompilerOutput = TRUE' to see C++ compilation details.
conf <- configureMCMC(Rmodel,monitors=params, thin=1, useConjugacy = TRUE)
## ===== Monitors =====
## thin = 1: D, id.prob, lam0.mark, lam0.resight, N, psi, sigma
## ===== Samplers =====
## slice sampler (21560)
##   - y.full[]  (21560 elements)
## RW sampler (1504)
##   - lam0.resight
##   - lam0.mark
##   - sigma
##   - id.prob
##   - s[]  (1500 elements)
## conjugate sampler (1)
##   - psi
## binary sampler (750)
##   - z[]  (750 elements)
#### --- Usamos ahora el MH  --- ####
conf$removeSampler("y.full")
for(j in 1:J.resight){
  for(k in 1:K){
    conf$addSampler(target = paste(paste("y.full[1:",n.marked,",",j,",",k,"]"), sep=""),
                    type = 'IDSampler',
                    control = list(nnid = nnid[j,k], j=j,k=k,n.marked=n.marked),
                    silent = TRUE)
  }
}

# Compilamos con el nuevo muestreador
Rmcmc <- buildMCMC(conf)


Cmcmc <- compileNimble(Rmcmc, project = Rmodel)
## Compiling
##   [Note] This may take a minute.
##   [Note] Use 'showCompilerOutput = TRUE' to see C++ compilation details.
nb = 1000
ni = 5000 + nb
nc = 3

# No ejecutamos el modelo en clase por el tiempo de ejecución
# outNim <- runMCMC(Cmcmc, niter = ni , nburnin = nb , nchains = nc,
#                   setSeed = FALSE, progressBar = TRUE, 
#                   samplesAsCodaMCMC = TRUE)

load("outNimZorro.RData")

Obtenemos resultados e inspeccionamos la convergencia

Podemos ver que aunque no se disponía de ningún identificado, se consigue una estima muy precisa (CV < 18%)

options(scipen=999)
summary(outNim)
## 
## Iterations = 1:5000
## Thinning interval = 1 
## Number of chains = 3 
## Sample size per chain = 5000 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##                    Mean        SD   Naive SE Time-series SE
## psi            0.294613  0.052054 0.00042502      0.0029844
## lam0.resight   0.061985  0.010039 0.00008197      0.0006391
## lam0.mark      0.009913  0.003257 0.00002659      0.0000746
## sigma          0.779872  0.015694 0.00012814      0.0003788
## N            220.487667 37.091549 0.30285123      2.2205418
## D              0.947353  0.159368 0.00130124      0.0095408
## id.prob        0.082871  0.076694 0.00062621      0.0020952
## 
## 2. Quantiles for each variable:
## 
##                    2.5%        25%       50%       75%     97.5%
## psi            0.207394   0.257762   0.28960   0.32588   0.41035
## lam0.resight   0.043910   0.054992   0.06145   0.06856   0.08252
## lam0.mark      0.004754   0.007607   0.00948   0.01176   0.01739
## sigma          0.749361   0.769416   0.77959   0.79003   0.81133
## N            159.000000 194.000000 217.00000 243.00000 304.00000
## D              0.683163   0.833545   0.93237   1.04408   1.30617
## id.prob        0.002365   0.025822   0.06061   0.11769   0.28690
xyplot(outNim)

Referencias

De Valpine, P., Turek, D., Paciorek, C. J., Anderson-Bergman, D., Lang, T., & Bodik, R. (2017). Programming with models: writing statistical algorithms for general model structures with NIMBLE. Journal of Computational and Graphical Statistics, 26, 403–413.DOI:10.1080/10618600.2016.1172487

Jiménez, J., Chandler, R. B., Tobajas, J., Descalzo, E., Mateo, R., & Ferreras, P. (2019). Generalized spatial mark–resight models with incomplete identification: An application to red fox density estimates. Ecology and Evolution, 9(8), 4739–4748. https://doi.org/10.1002/ece3.5077

R Core Team. (2022). R: A language and environment for statistical computing. R Foundation for Statistical Computing. Vienna, Austria. Retrieved from https://www.r-project.org/.